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Abstract 

We develop a high-order kinetic scheme for entropy-based moment models of a one-dimensional linear kinetic 
equation in slab geometry. High-order spatial reconstructions are achieved using the weighted essentially 
non-oscillatory (WENO) method, and for time integration we use multi-step Runge-Kutta methods which 
are strong stability preserving and whose stages and steps can be written as convex combinations of forward 
Euler steps. We show that the moment vectors stay in the realizable set using these time integrators along 
with a maximum principle-based kinetic-level limiter, which simultaneously dampens spurious oscillations in 
the numerical solutions. We present numerical results both on a manufactured solution, where we perform 
convergence tests showing our scheme converges of the expected order up to the numerical noise from 
the numerical optimization, as well as on two standard benchmark problems, where we show some of the 
advantages of high-order solutions and the role of the key parameter in the limiter. 
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1. Introduction 

In recent years many approaches have been considered for the solution of time-dependent linear kinetic 
transport equations, which arise for example in electron radiation therapy or radiative heat transfer problems. 
Many of the most popular methods are moment methods, also known as moment closures because they are 
distinguished by how they close the truncated system of exact moment equations. Moments are defined 
through angular averages against basis functions to produce spectral approximations in the angle variable. 
A typical family of moment models are the so-called P^v-methods [16, 25] which are pure spectral methods. 
However, many high-order moment methods, including Pjv, do not take into account that the original 
kinetic density to be approximated must be nonnegative. The moment vectors produced by such models are 
therefore often not realizable, that is, there is no associated nonnegative kinetic density consistent with the 
moment vector, and thus the solutions can contain obviously non-physical artifacts such as negative local 
particle densities [6]. 

The family of minimum-entropy models, colloquially known as Mjv models or entropy-based moment do- 
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sures, solve this problem (for certain physically relevant entropies) by specifying the closure using a non¬ 
negative density reconstructed from the moments. The Mjv models are the only models which additionally 
are hyperbolic and dissipate entropy [24] . The cost of all these properties is that the reconstruction of this 
density involves solving an optimization problem at every point on the space-time mesh. These reconstruc¬ 
tions, however, can be parallelized, and so the recent emphasis on algorithms that can take advantage of 
massively parallel computing environments has led to renewed interest in the computation of Mjv solutions 
both for linear and nonlinear kinetic equations [2, 10, 15, 18, 21, 27]. Despite the parallelizability of the 
cost of the numerical optimization, the gain in efficiency that would come from a higher-order space-time 
discretization will still be necessary for a practical implementation. 

The key challenge for high-order methods for entropy-based moment closures is that the numerical solutions 
leave the set of realizable moments [28] , outside of which the defining optimization problem has no solution. 
Discontinuous-Galerkin methods can handle this problem using a realizability limiter directly on the moment 
vectors themselves [4, 28, 39], but at this level realizability conditions are in general quite complicated and 
also not well-understood for two- or three-dimensional problems for moment models of order higher than two. 
Realizability limiting for kinetic schemes, however, is much easier because at the level of the kinetic density, 
realizability corresponds simply to nonnegativity. Furthermore, this same limiter can be strengthened to 
also enforce a local maximum principle, thereby dampening artificial oscillations in numerical solutions. 

Thus in this work we derive a high-order (in space and time) kinetic scheme for Mjy models with moments 
of (in principle) arbitrary order. We start in Section 2 by introducing the linear kinetic equation we will con¬ 
sider, its entropy-based moment closure, and reviewing the concept of realizability. Continuing in Section 3 
we introduce the concept of a kinetic scheme for moment equations and then give our numerical techniques 
for the discretization of each of the independent variables: angle, space, and time. The issue of realizability 
preservation and the necessary limiters are discussed in Section 4, finishing the full description of our scheme. 
The results from our numerical simulations are presented in Section 5, including a convergence study using a 
manufactured solution and solutions for two benchmark problems. Finally we draw conclusions and discuss 
the next steps for future work in Section 6. 


2. A linear kinetic equation and moment closures 

We begin with the linear kinetic equation we will use to test our algorithm and a brief introduction to 
entropy-based moment closures which closely follows [4]. More background can be found for example in 
[18, 24, 25] and references therein. 

2.1. A linear kinetic equation 

We consider the following one-dimensional linear kinetic equation for the kinetic density ip = ip(t,x,^) > 0 
in slab geometry, for time t > 0, spatial coordinate x G X = (a:L,a:R) C K, and angle variable /i G [—1,1]: 

dtip + fj.dxip + (T^ip = asC{ip) + S, (2.1) 

where da and Us are the nonnegative absorption and scattering interaction coefficients, and S a source. The 
operator C is a collision operator, which in this paper we assume to be linear and have the form 

C(V’)=y T{fj.,fi')ip{t,x,fj.') dn'- J T{fi',fi)ip{t,x,fi) dfj.'. (2.2) 

We assume that the kernel T is strictly positive and normalized to T{ii', lAjdfj,' = 1. A typical example 
is isotropic scattering, where = 1/2. 
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Equation (2.1) is supplemented by initial and boundary conditions: 


1p{t,Xi,,p) = 

= V’L(i,M), 

t > 0, 

p>0, 

(2.3a) 

1p{t,XB.,p.) = 


t > 0, 

/r < 0, 

(2.3b) 

ilj{0,x,p.) = 

= 1pt=o{x,fl) , 

X e (xL,a;R), 

fi € [—1,1] , 

(2.3c) 


where V’Li V'Rj ipt=o are given. 


2.2. Moment equations and entropy-based closures 

Moment equations are an angular discretization for (2.1), where the moments themselves are defined by 
angular averages against a set of basis functions. We use the following notation for angular integrals: 

{cj)) = j (j){p)dp 

for any integrable function cj) = (/<(/r); and therefore if we collect the basis functions into a vector m = 
m(^) = (mo(p.), mi(/i),... ,mAr(/i))^, the moments of a kinetic density 4> are given by u = (m^). In this 
paper we consider the monomial moments mi{pL) = ijd, though all results can be extended to other bases, 
including, for example, partial [11, 12] or mixed moments [13, 32]. 

The closed system of moment equations is a system of partial differential equations of the form 

dtu + 5a;f(u) + (TaU = crsr(u) + (mS ), (2.4) 

where the moment vector u(t, x) approximates {rxiijj) for the kinetic density ip satisfying (2.1). In an entropy- 
based closure (commonly referred to as the Mat model or the Levermore closure after he exposed their general 
structure in [24]), the functions f and r have the form 

f(u) := (^fj.niipjj and r(u) := (^inC{ipu)'j . 

Here V'u is an ansatz density reconstructed from the moments u by solving the constrained optimization 
problem: 

V'u = argmin{(? 7 (((i)) : (m^) = u} , (2.5) 

4 > 

where the kinetic entropy density r] is strictly convex and the minimum is simply taken over functions 
(j) = 4 >{pl) such that {ri{(j))) and {mcp) are well defined. This problem is typically solved through its strictly 
convex, unconstrained, finite-dimensional dual, 

q:(u) := argmin (77*(m^a)) — u^a, (2-6) 

a,gR«+i 

where ry* is the Legendre dual of rj. The first-order necessary conditions for q:(u) show that the solution to 
(2.5) has the form 

= r]'^ (m'^Q:(u)) (2.7) 

where is the derivative of 77 *. 

The kinetic entropy density rj can be chosen according to the physics being modeled. 

While in a linear setting such as ours, indeed any convex entropy r] is dissipated by (2.1). As in [18] we 
focus on the Maxwell-Boltzmann entropy. 


r]{z) = z\og{z) - z, 
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because not only is it physically relevant for a wide variety of problems, but in particular it gives a positive 
ansatz ip-a, since r 7 *( 2 /) = ri'^{y) = exp{y), and thus 

'0U = exp (m^Q;(u)) . (2.8) 

Another standard choice for the entropy is ri{z) = which yields the well-known Pjv equations [18, 
25]. The Pat equations are linear, and since = y, the optimization problem (2.6) can be solved by 

hand. However the resulting ansatz is simply a linear combination of the basis polynomials and thus is not 
necessarily nonnegative. ^ 

The incorporation of the boundary conditions (2.3) is neither obvious nor trivial and is still an open problem 
[22, 23, 29, 35]. This is not a focus of our work here, and therefore we only use a simple approach from 
previous work on entropy-based moment closures [2, 18], which we discuss below in Section 3.2. 


2.3. Moment realizability 

Since the underlying kinetic density we are trying to approximate is nonnegative, a moment vector only 
makes sense physically if it can be associated with a nonnegative density. In this case the moment vector is 
called realizable. Additionally, since the entropy ansatz has the form (2.8), the optimization problem (2.5) 
only has a solution if the moment vector lies in the ansatz moment space 

A := {(mexp (m^a)) : a G . 

In our case, where the domain of angular integration is bounded, the ansatz moment space A is exactly 
equal to the set of realizable moment vectors [19]. Therefore we can focus simply on realizable moments: 

Definition 2.1. The realizable set TZ^ is 

T^m = • 30(^) > 0, {(j)) > 0, such that u = (mi)))} . 

Any (j) such that u = {nKp) is called a representing density. 


The realizable set is a convex cone. In the monomial basis, a moment vector is realizable if and only if its 
corresponding Hankel matrices are positive definite [9, 33]. 

In general, angular integrals cannot be computed analytically. We define a quadrature for functions (j) : 
[—1,1] —)■ M by nodes and weights such that 

Nq 

Wi4>{y,i) « {(t>) 

i=l 


Then the numerically realizable set is [3] 


Nq 


TV^ = 


3fi > 0 S.t. U = ( C 


i=l 


Indeed, when replacing the integrals in the optimization problem (2.5) with quadrature, a minimizer can 
only exist when u G 7?.®. Below we often abuse notation and write (</>} when in implementation we mean 
its approximation by quadrature. We also liberally use the term realizable either to mean realizability with 
respect to TZ^ or with respect to TZ^, where the specific meaning depends on whether exact integrals or 
those approximated by quadrature are meant in the context. 


^ For this reason, while the rest of the scheme we give below can be applied to the Pat equations, the positivity-preserving 
techniques we use (see Section 4) would be unnecessary. However, the scheme should apply equally well for any entropy-based 
closure with a positive ansatz. 
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3. A high-order kinetic scheme 


A kinetic scheme for (2.4) can be thought of as first defining a spatial discretization for the underlying kinetic 
equation (2.1) and subsequently performing the angular discretization with the moment closure. We divide 
the spatial domain (xljSJr) into a (for simplicity) uniform grid of J cells Ij = {xj-ij 2 ,Xj+i/ 2 ), where the 
cell edges are given by Xj±i /2 = Xj ± Aa:/2 for cell centers Xj = xr + (j — l/2)Ax, and Ax = (xr — xl)/J- 
For (2.1) we define a finite-volume scheme for the cell means 

'i/’(t, X, /r)dx, (3-1) 

which with the Godunov (or ‘upwind’) numerical flux gives: 

dtipj + max(/i, 0) 

+ min(/i, aaipj = ^asC{'4))- + Sj 

where V’^ 1/2 terms denote the values of the approximate solution at the cell edges 

Xj±i /2 from the left and right, respectively, and we generally use the bar with subsequent subscript j, i.e. 
^j, to indicate a cell average over the j-th cell as in (3.1). 

To obtain a high-order scheme in space one only has to give a high-order reconstruction of the point-values 
of ip, the distribution underlying the cell means, not only at the edge values for the flux terms, but also 
throughout the cells when Ca or dg depends on x. We use the popular weighted essentially non-oscillatory 
(WENO) reconstruction method [34], which gives a polynomial reconstruction of ip from the cell averages 
of the j-th cell and its neighbors. 

Now we perform the moment closure by replacing '0 with the entropy ansatz ip in (2.7), multiplying through 
by the angular basis functions, and integrating out the angle: 



dtu.j + ^/rm' 


V’,+1/2 - 


f-1/2 


Ax 


/rm 


V ’j- + l /2 - V’j-1/2 ^ 

Ax , 


CTaUj = 2^sr(u)^- 


(3.2) 


where 





and 


(</>)- 




/-I 


and s = (mS'). Here 0^i/2 denote the evaluations at the cell edges of the WENO reconstructions made 
using the entropy ansatze of the neighboring cells evaluated from the right, and respectively for 0 ’~j_ 2/2 
evaluated from the left. The first step in the scheme, then, is to compute the ansatz for each cell. 


3.1. Numerical optimization for angular reconstruction 

In order to compute ipj at the cell means, we first compute the multipliers Q:(uj) by solving the dual problem 
(2.6). The gradient and Hessian of the objective function are 

g(Q:) = (mexp(m’^Q:)) — Uj and H(a) = (mm^ exp(m’^a)) . (3-3) 

We use the numerical optimization techniques proposed in [3]. 
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We assume that the moment vector has been scaled so that its zeroth component is one. The optimizer 
stops at the first iterate a which satishes 


l|g(a)ll < 1 , IIJ II , =: t', and (3.4a) 

1 + \m\\ +T 

1-e < exp(-||d(a)|li - | log(uo(a))|) (3.4b) 


where || • || is the Euclidean and || -Hi the 1—norm in r and e user-specified tolerances, and uo{a.) = 

(exp(m^Q:) is the zero-th order moment associated with the multiplier vector a. This criterion is similar to 
the one in [3] but modified for the following reasons. 

Unlike the algorithm there, we modify the zero-th component of the final multipliers so that the zero-th 
moment (which gives the local density) is exactly matched. That is, while the optimizer stops at the 
first a which satisfies (3.4), the multiplier vector it returns (to define 'ipj for the kinetic reconstructions in 
Section 3.2) is 

a := (ao - log(Mo(Q:)),Q!i,... ,Q;Ar) . (3.5) 

Since mo = 1, this gives (mexp(m^a)) = (mexp(m'^a))/Mo(a), which ensures (exp(m^a)) = 1. The 
form of t' on the right-hand side of (3.4a) ensures that ||g(a)|| is bounded by r: 


||g(a)|| = II (mexp(m^a)) - Uj 
1 


(3.6) 


< 


< 


uo{a) 
1 

1 — t' 


(mexp(m'^a)) — Uj|| -f 


uo{a) 


1 — r' 


U, = T, 


where we have used that (3.4a) implies |uo(q!) — 1| < t' and consequently 1 /uo(q: )<l/(l-r'). 

The second condition (3.4b) in the stopping criterion enforces an approximate lower bound on the ratio^ 
where = exp(m^Q:(uj)) is the entropy ansatz associated with an exact solution of (2.8) 

and ipj = exp(m’^cS) is the ansatz associated with the Lagrange multiplier vector which satisfies (3.4). While 
of course we do not know the exact entropy ansatz, we can approximate the difference between the exact 
solution Q;(uj) and another multiplier vector ct from the iterations of the optimizer by the Newton direction 

d(a) = -H(Q:)"^g(a). 


This leads to the approximate bound 

■0 ■ 

^ = exp(m^(Q:(uj) — a)) = exp(m^(d(uj) — a. + a. — a)) 

« exp(m^(d(a) + a. — «)) 

> exp(-||m||oo(||d(Q;)||i -k I log(uo(a))|)), 

where |jm||oo = max^,^ \'mi{n)\ = 1. 

Finally, and exactly as in [3], we use a isotropic-regularization technique to return multipliers for nearby 
moments when the optimizer fails (for example, by reaching a maximum number of iterations or being unable 
to solve for the Newton direction). Isotropically regularized moments are defined by the convex combination 

v(u,r) := (1 - r)u-k ruoUiso, (3.7) 


^This is also similar to what was done in [3], though there the bound was more naturally given as an upper bound of the 
inverse of the ratio we use. 
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where Uiso = 5 (m) is the moment vector of the normalized isotropic density = 1/2. Then the optimizer 
moves through a sequence of values 0 ,ri,r 2 ,... ,rM, advancing in this sequence only if the optimizer fails 
to converge for v(u, r) after k-ceg iterations for the current value of r. It is assumed that tm is chosen large 
enough that the optimizer will always converge for v(u, tm) for any realizable u. 

3.2. Spatial WENO reconstruction 

We use the standard WENO reconstruction method given for example in [34, 36]. For the unfamiliar 
reader, in this section we briefly introduce the method, while a more detailed documentation and demo 
implementations of the reconstruction procedures we implemented can be found on our webpage [ 1 ]. 

Here, as in the previous section, we use ipj to indicate the entropy ansatz associated with the Lagrange 
multipliers returned by the optimizer to approximate the true multipliers Q;(uj). Time dependence is again 
suppressed for clarity of exposition. 

At a: = Xj-i/ 2 , the cell edge between the (j — l)-th and j-th cells, for each p, we evaluate a weighted 
combination of polynomials of degree fc — 1 , Pjm{.‘,ti) G Pfe-i (the space of polynomials up to degree k — 1 ), 
m = 0 ,l,...,fc, each solving the interpolation problem 

Pjmix,p)dx = tlj,,ip), £€ {j -k + m,...,j + m-l}. (3.8) 

The WENO method then gives weights m form the weighted averages 

k 

m—0 

and finally we approximate the values at the cell interfaces by 

^ 7 _i/ 2 (m) -P 7 _i/ 2 (a^i-i/ 2 ,M), and ^/’+_i/ 2 (/r) ^ p'^_^/^{xj_i/ 2 , p). 

The weights wj 7 i /2 m non-linear functions of the cell-averages and reflect the smoothness of each poly¬ 
nomial Pjm- They are computed such that for smooth data the approximation order at the cell edge is 
maximized. This gives an order 2fc — 1 approximation at the cell edge, while the overall order in the interior 
of the cell is k. 

When at least one of the interaction coefficients da or (Tg is spatially dependent, we must also specify the 
reconstruction inside each cell to compute, for example, the daU^ term in (3.2). Here we must make a choice, 
because both P^_i /2 ^*7+1/2 order k reconstructions of the density ip in the j-th cell.^ We denote 

this polynomial 'ipj{x,p) and choose it to be 

{ P~+i/2i^^p) itp>0, 

i (p7+i/2(^>^)+^^- 1 / 2 ( 2 ^’^)) if 1 ^ = 0. (3-9) 

p 7 -i/ 2 (^>I^) ifAi< 0 . 

This particular reconstruction allows us to derive a realizability-preserving time step in Theorem 4.1. 

Finally, it remains to incorporate boundary conditions. We define ‘ghost cells’ at the cell indices j G 
{I — k,... ,0, J + 1,... ,J + k}, namely those indices which are used in (3.8) but have not yet been defined. 


® Some reconstruction methods, such as subcell WENO [7] or minmod [36], give only one polynomial reconstruction inside 
each cell, and so for these methods such a choice would be unnecessary. 
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We assume that we can smoothly extend V’l and V'R in M to [—1,1]- We then use the simplest possible 
approach^ and set 


if j e {l-A:, ...,0} 


(3.10a) 


We note, however, that the validity of this approach is not entirely noncontroversial, but the question of 
appropriate boundary conditions for moment models is an open problem [22, 23, 29, 35] which we do not 
explore here. 

In case of periodic boundary conditions, we use the data from the physical cells j S {J — k + 1,..., J} from 
the right side of the domain to fill the ghost cells j G {1 — k,..., 0}. The ghost cells on the right side of the 
domain analogously take the values from the left side of the physical domain. 

3.3. High-order time integration 

If we collect the approximate cell means from each spatial cell into one long vector Uh{t) := (u^ (!)> (i)j • ■ • > u 

then equation (3.2) can be written as 

dt^h = Lh{uh) 

Since entropy-based moment closures are only defined on the realizable set, it is important to choose a time 
integrator for which we can prove that realizability is preserved. Therefore we follow [2, 39] and use a strong 
stability-preserving (SSP) method whose stages and steps are convex combinations of forward Euler steps. 
Since the realizable set is convex, the analysis of a forward Euler step then suffices to prove realizability 
preservation of the high-order method. 

When possible we use a SSP Runge-Kutta (SSP-RK) method, but such methods only exist up to order four 
[17, 31]. Eor higher orders we use the so-called two-step Runge-Kutta (TSRK) SSP methods [20] as well as 
their generalizations, the multi-step Runge-Kutta (MSRK) SSP methods [5]. They combine Runge-Kutta 
schemes with positive weights and high-order multistep-methods to achieve a total order which is higher 
than four while preserving the important SSP property. If we let indicate the collection of the numerical 
approximations to the cell averages of the solution at the n-th time instant = nAt, an s-stage TSRK 
method in the low-storage implementation has the following form [ 20 ]: 


Vi — ^ -|- ^1 

ur' = cur^ + (^i- 

where the coefficients define the scheme. In this work we only consider explicit schemes, 

where gim = 0 for £ > m. The positive coefficient p is called the radius of absolute monotonicity and 
indicates how much we can scale our time step At while fulfilling the CEL condition for forward Euler steps 
(see Section 4.1 below). 

Such schemes provide reasonably good effective CEL numbers, which are the ratios p/s for each method. 

If the forward-Euler method is stable under the timestep restriction At < Atp, the high order scheme with 
radius of absolute monotonicity p is stable under the time-step restriction At < pAtp. A larger effective 


At 


d^ ^ ^ Qim J T 'y ^ Qim j Vm Lh^Vm) j 5 0 A ^ ^ 


m—0 


m—0 


At 


C y Vm j -|- y ] r]m I Um + LhiUm) j > 
m—0 j m—0 \ ^ J 


a smooth setting this might reduce the order of accuracy at the boundary to one. 




TS 

[IK 

TS 

=IK 


3SPRK 
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th 


Figure 1 - One possible startup procedure for SSP TSRK schemes. The first step from to to ti is subdivided 
into substeps (here there are three substeps of sizes At/4, At/4, and At/2). A one-step SSP Runge-Kutta 
scheme is used for the first substep, and subsequent substeps are taken with the TSRK scheme itself, doubling 
the step sizes until reaching ti. Illustration taken from [20]. 


Order Method(TO, k, s) Effective CEL 


1 

SSPRK(1,1,1) 

1.0 

2 

SSPRK(1,2,20) 

0.95 

3 

SSPRK(1,3,16) 

0.75 

4 

SSPRK(1,4,10) 

0.6 

5 

TSRK(2,5,8) 

0.4474 

6 

TSRK(2,6,12) 

0.3653 

7 

MSRK(5,7,12) 

0.3089 


Table 1 - Methods used here and their effective CFL. The nomenclature METHOD(m, k, s) denotes a method 
with m steps, order k, and s stages. 


CEL indicates that, in order to reach a given final time, the operator Lfi will need to be evaluated fewer 
times. The explicit Euler method, which serves as the reference for efficiency in this context, has an effective 
CEL condition of 1 (s = p = 1). 

The effective CFL numbers of the integration schemes used here are given in Table 1. 

Unfortunately, multistep methods are not self-starting, so they need a predictor for the first time-step. 
Since we can only use convex combinations of forward Euler steps for all time steps in order to prove that 
realizability is maintained we must use a lower-order method. We use the strategy given in [20]: First, we 
predict with a smaller step size At* = At/2'J, for an integer q> 1, with the ten-stage, fourth-order explicit 
SSPRK(1,4,10) method. Then we use the corresponding TSRK method and double the step-size after every 
iteration until we reach t = At. This procedure is shown in Figure 1 for q = 2. 

For the five-step method MSRK(5, 7,12), which we use for seventh-order simulations, we have to initialize 
four steps. For these initialization steps we use the two-step method TSRK(2, 7,12), whose initial step we pre¬ 
dict using the same method given above. However the radius of absolute monotonicity p for TSRK(2, 7,12) 
is approximately 2.7659, while for MSRK(5, 7,12), the radius of absolute monotonicity p is approximately 
3.0886. This means the time steps we take after initialization will be longer than those we can take with the 
TSRK method during initializating without violating the realizability-preserving CFL condition. Therefore 
we stop increasing the step size in the initialization routine when we reach At/2 (as opposed to At), and 
then continue with the TSRK initialization steps of size At/2 until we have computed Ufi(t 4 ) = Uft,(4At). 
At this point we have all the previous steps we need to compute u^(t 5 ) using the five-stage MSRK(5, 7,12) 
method. 


4. Realizability preservation and limiting 

The strategy to prove that the moments Uj remain realizable at each time step and inner stage of the time 
integrator begins by proving that the kinetic density reconstructed for the cell means remains nonnegative 
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after a forward Euler step with a certain time-step restriction. Since we use time integrators which are a 
convex combinations of Euler steps (see Section 3.3) this immediately gives nonnegativity for all time steps 
and internal stages of time integration. This proof, however, requires the assumption that the point-wise 
values of the current polynomial reconstruction are nonnegative at every spatial and angular quadrature 
point. One therefore introduces a limiter to enforce this nonnegativity. 


4-1- Realizability preservation of the cell means 

We follow along the lines of the main proof in [2] and will provide weaker conditions which follow [37, 39]. 

A spatial quadrature rule plays a crucial role here. We use Gauss-Lobatto rules which are exact for polyno¬ 
mials of degree 2Q — 3, where Q is the number of quadrature nodes. These rules are characterized by nodes 
Pi and weights Wi for i G { 1 , 2,..., Q} on the reference interval [—\^ 5 ]: 

|■l/2 Q 

/ f{x) dx (4.1) 

^-1/2 

We let Xji := Xj + Axpi denote the quadrature nodes shifted and scaled for cell Ij, and note that since the 
Gauss-Lobatto rules include the endpoints, we have xji = Xj- 1/2 and xjq = Xj+i/ 2 - 

We also use the property that, since the collision kernel T in (2.2) is nonnegative, there exists a realizable 
moment vector such that 

(mC(V’)) = uc - u. (4.2) 

Theorem 4.1 (Main theorem). Assume that 


(i) for all cells j G {1, 2,..., J} we have 0 < S{tn,x)\i., cra{tn,x)\i., and crs{tn,x)\i. are in Fks-i{Ij); 

(ii) the cell means u" at time step t^ are realizable; 

(Hi) in each cell the ratio between the exact entropy ansatz and its approximation from the optimizer satisfies 


A Ad) ~ 


\ — e 


for all fj, G [— 1 , 1 ] in the angular quadrature set; and 


(4.3) 


(iv) 


the point-wise values of the polynomial reconstructions ipj{x,pL) G P/c-i at the quadrature nodes of the 
Q-point Gauss-Lobatto on each cell Ij are nonnegative, for Q = |"(fc -I- fcs -I- l)/2] 


Let 


^t,max •— max max cy^{t,^,Xjf) -t- crXjf). 


(4.4) 


Then under the CFL condition 


At < (1 — e) min 


Axwr 


<^t,max 1 -I- Axwqat ,r 


the cell means Uj after one forward Euler step are realizable. 


(4.5) 


^Where f-] is the ceiling function, that is, it returns smallest integer bigger than or equal to its argument. Since the Gauss- 
Lobatto rule is exact for polynomials of degree 2Q — 3 this choice guarantees to exactly integrate the occurring polynomials of 
degree {k kg — 2). 
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Proof. For simplicity we will neglect the time-index for quantities from t and use it only for time-level 


t = tn+i- An Euler step is given by 

^ " "^^”- 1 / 2 ) )+ + ('^/+i/2 - ^U/ 2 ) )_) (4-6a) 

At -I- crsr(u)j- -I- Sj^ (4.6b) 

and consequently we have = {^ 4 >j) for (fj given by 

^ ^ (max(/r, 0) (V^’+i/a “ ^J- 1 / 2 ) + 0) (V^Z+i/a “ '^t- 1 / 2 ) ) 

-I- At {-ut'ipj + (JsilJCj + Sj) (4.7b) 


where the total interaction coefficient is defined as cxt := da + CTs, and tfc ^ 0 is the entropy ansatz 
corresponding to the realizable part uc of the collision operator (see (4.2)). Note that (fj = depends 

on /X. 

Let us first consider the case /x > 0. Stripping away positive terms and using /i < 1 gives 

- Af _ _ 

- ^ 3 + 1 / 2 - (4.8) 

Next, we want to use our polynomial reconstruction (3.9) and the Q-point Gauss-Lobatto quadrature to 
relate the exact cell mean '0^ to the reconstruction’s values at the cell edge 0 ^ 01/2 average Uti^y 

But since only the approximations from the optimization are used to define the reconstruction, we must 
first multiply and divide ipj by tpy Note that this is possible since ipj > 0. Then we use the assumed bound 
(4.3), note carefully that for /x > 0 indeed ~ '^ji^jQ))^ apply L°° bounds on the total interaction 

coefficient. This gives 


( 4 . 8 ) ^ ^ At ^ 

<t>j > wYjiXji) - - At^ Wiat{Xji)%fj{xy) 

i=l 


Aa; 


i=l 


( 4 . 3 ) 


Af 


- - At'^WiCrt{Xji)lpj{Xji) 

( 4 . 4 ) [ At \ 

> ^ xii (1 - £ - Af(Tt,max) 'f’ji + 1 ( 1 - SJWQ - ^ - Atwgat ,max V’aQ, 


where in the last line we have introduced the notation ^pji := ifj^Xji). One can see that (4.5) ensures 
nonnegativity of both terms in the final expression. Recalling that xui = xig, the cases /x < 0 and /x = 0 
follow analogously, and together we have that (fj > 0 at each fj,, which shows that is realizable. □ 

Remark 4.2. In practice, (4.3) is only approximately enforced by (3.4). However, we have never had 
problems losing realizability in our numerical simulations, including those that approach the boundary of 
realizability. 

Remark 4.3. We have assumed that the reconstruction (3.9) is always used to compute uYj- However, 
when (Ta and are constant in each cell, one simply has cr^j — therefore the reconstruction 

is not needed. The proof can then be completed with an upper bound on ifj/ipj, which is similarly easy to 
enforce approximately, and one ends up with a slightly less restrictive CFL condition. However, we did not 
use this in our implementation. 
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The time-step restriction (4.5) can then be scaled by p for the corresponding time-integration scheme to 
give realizability of every stage and step in the scheme. 

In practice, to achieve an order k method for sources S or interaction coefficients Ua or Us which are not 
piecewise degree k — 1 polynomials, one would approximate them using the same spatial reconstruction 
techniques that we use for the density to achieve an order k approximation of the corresponding terms. 
Thus in Theorem 4.1, one would not use a value of fcg larger than k. 


4-2. Limiting 

The first role of the limiter, then, is first to ensure that the point-wise values of the polynomial reconstructions 
are nonnegative at the spatial and angular quadrature points. However, as we will see in Section 5.3 
below, the numerical solutions using a limiter which only ensures nonnegativity can still contain spurious 
oscillations. Therefore we extend the same limiter to enforce local maximum principles as well, thereby 
much more effectively dampening such oscillations. 


4-2.1. Positivity-preserving limiter 


To preserve nonnegativity we can simply apply a linear scaling limiter. The limited spatial reconstruction 
is defined as tjjjix) := (1 — + Otpj{x)\ notice that 0 = 1 corresponds to no limiting. For each quadrature 

point (both in space and angle, though here we suppress the angular argument) we compute 


_ 

i’j - i’ji 

1 


if < 0, 


else. 


(4.9) 


Then in each cell we set 

^ = ^3 ■= . min J 

(where one should keep in mind that 6 j still depends on p). One immediately sees that this limiter ensures 
the positivity, preserves the cell means ipj, and following arguments from [37, 38], does not destroy accuracy 
of the scheme if tjjj > 0- 

However, it has been remarked in [41] that in some pathological situations this limiter may reduce the 
accuracy to second order. See also [4] where a similar observation has been made for a realizability limiter 
in a discontinuous-Galerkin scheme. 


4-2.2. Maximum principle-satisfying limiter 

The limiter we introduce here is a slightly modified version of the maximum-principle limiter from [40]. 

Since we know a priori that ip satisfies a strict maximum principle m < ip{x.,p) < M for all x and /x, 
a natural strategy to dampen artificial oscillations in numerical solutions is to enforce a local maximum 
principle. Specifically, we would like the polynomial reconstruction ipj{x) to be bounded by the data of 
those cells which influence it. The corresponding index set of influential nodes is = {j ~ ..., j -\- k} 

(cf. (3.8)), so the local maximum principle we would like to enforce is 

min ipeip) <ipj{x,p) < max 

However this tends to flatten smooth extrema, so, inspired by the modified minmod function in [8] , we relax 
the strict maximum principle by setting the maximum principle bounds Mj and ruj locally as 

/ \ _ / ^x \ _ 

Mj := I 1 -\-c - I max ipAu) and m, := 1 — c- I min ipAu), 

^ \ 2 y ^ V 2 y ^ 
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where c is a local bound on the relative derivative of tp, i.e. max^^a; \dx'tp{x,^)/il}{x^pL)\^ and the maximum 
and minimum in /i are taken over the angular quadrature nodes.® Therefore the maximum principle that 
we will actually enforce is 

xnj < ipj {xji , /i) < Mj 

for all spatial quadrature points Xji S Ij and all angular quadrature points. To enforce this maximum 
principle, for each spatial quadrature point we set 


Mj - iPj 


Qji :— < 


P’21 - 
rrij — Ipj 

P’ji - P’j 

1 


^ \iipji>Mj, 


if tpji < rrij, 
otherwise, 


(4.10) 


and finally for each cell we choose 9j := mmi{9ji}. 

As with the positivity-preserving limiter, it can be shown that this limiter does not destroy accuracy 


5. Numerical results 

In this section we present results to confirm that our scheme converges with the expected order, to show the 
effect of various parameters in the scheme, and to highlight some of the features of high-order solutions. 

Except where otherwise noted, we used the following parameter values: 

r = 10“® Optimization gradient tolerance, 

e = 0.01 Optimization tolerance on 1 — ipj/ipj, 

{ri ,..., xm} = {10“®, 10“®, 10“'*} Outer regularization loop in optimizer, 

kr = 50 Number of optimization iterations before 

advancing outer regularization loop, 
riQ = 40 Number of angular quadrature nodes, 

c = 1 Bound on \dxip/ip\ in maximum-principle limiter. 

For the angular quadrature we used (NQ/2)-point Gauss-Lobatto rules over both /r G [—1, 0] and /i G [0,1]. 

For the value q determining the number of initialization steps for the multi-step time integrators we used 
two for fifth-order simulations and three for sixth- and seventh-order simulations. 

The time step is chosen to fulfill (4.5) (replacing At by At/p for the appropriate time integrator) with 
equality. 

In both benchmark problems we use isotropic scattering, C{ip) = \{ip) — ip. 

For the first stage of the first time step, the initial multipliers for the optimizer are those associated with the 
normalized isotropic distribution. For the following stages and time steps, the initial multipliers are set to 
those from the previous stage or step at the same spatial cell. However, if in these later stages the optimizer 
cannot converge before initializing the regularization loop, we first switch back to the multipliers associated 
with the normalized isotropic distribution and restart the optimizer (this time allowing regularization if 
necessary). In our experience, the isotropic multipliers are the safest choice for the initial condition, and 
therefore this technique reduces the number of times the regularization must be used. 


® Note that we always want mj > 0, so we never take c > 2l/\x. 
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As in [ 2 ], if regularization has to be applied to u^, we replace it (in (4.6)) by the regularized moments 
v(uj,r) (3.7) for which the optimizer converged. Then Theorem 4.1 can be applied to ensure that the next 
iterate will be realizable as well. 

In all figures below we plot the zeroth-order component of the reconstruction Uj(a;) = fi)) for x ^ Ij, 

see (3.9). 

5.1. M]\[ manufactured solution 

In general analytical solutions for minimum-entropy models are not known. Therefore, to test the conver¬ 
gence and efficiency of our scheme, we use the method of manufactured solutions, and we follow the target 
solution given in [4] but add a spatially and temporally dependent absorption interaction coefficient. The 
solution is defined on the spatial domain X = (—tt, tt) with periodic boundary conditions. 

A kinetic density in the form of the entropy ansatz is given by 

(j){t, X, = exp(ao(i, x) + ai{t, x)^), (5.1a) 

ao{t,x) = — K — sin{x — t) — a, (5.1b) 

ai{t,x) =K + sm{x — t). (5.1c) 

A source term is defined by applying the transport operator to (p: 

S{t,x,^j.) := dt<p{t,x,fj,) + nd^(j){t,x,fj,) + as,{t,x)(j){t,x,fj,), 

where 

cra(t, x) := 4 — 4cos(x — t)). 

Thus by inserting this S into (2.1) (and taking Ug = 0) we have that p is, by construction, a solution of 

( 2 . 1 ). 

A straightforward computation shows that S > 0 (for any a or K), which means that Theorem 4.1 will 
apply to the resulting moment system. Furthermore we take 

so that the maximum value of (p) for (t, x) G [0, tf] x X is one. As K is increased, p converges to a Dirac 
delta at /r = 1 . 

The moment vector w = {mp) is then a solution of (2.4) for Mw models with N > 1. 

We used the final time tf = 7 r /5 and chose K = 4, for which wi/wq G [0.67,0.8] (recall that \wi/wo\ < 1 
is necessary for realizability). We are using a fairly low value of K because, in order to show convergence 
for our scheme with the highest-order (fc = 7), we need to be able to have a tighter control on the errors 
from the numerical optimization. When K is higher, these errors are larger and drown out the convergence 
in space and time. In the following, we used the M 3 model so that our results included the effects of the 
numerical optimization. 

We compute errors in the zero-th moment of the solution, which we denote wo{t,x) = {p(t,x, •)). Then 
and errors for uo^h{t,x) (that is, the zero-th component of a numerical solution u/i) are defined as 

/ \wo{tf,x) - uo^h{tf,x)\dx and = ma.x\wo{tf,x) - uo^h{tf,x)\ 

Jx 

respectively. We approximate uo^h{t{,x) using the same reconstruction technique as for the scheme for the 
underlying kinetic density (3.9) and integrate with respect to /r. Then we approximate the integral in 
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k = 2 


CO 

II 1 


fc = 5 


II 1 


K 

V 

El 

y 

El 

y 

El 

y 

10 

6.532e-02 

— 

1.668e-02 

— 

4.816e-03 

— 

2.696e-03 

— 

20 

1.981e-02 

1.7 

9.931e-04 

4.1 

3.622e-05 

7.1 

5.130e-06 

9.0 

40 

3.823e-03 

2.4 

5.531e-05 

4.2 

2.517e-07 

7.2 

3.452e-09 

10.5 

80 

1.005e-03 

1.9 

6.808e-06 

3.0 

7.546e-09 

5.1 

5.090e-ll 

6.1 

160 

2.193e-04 

2.2 

9.778e-07 

2.8 

2.427e-10 

5.0 

4.049e-ll 

0.3 

320 

5.784e-05 

1.9 

1.317e-07 

2.9 

4.636e-ll 

2.4 

5.645e-ll 

-0.5 



V 

Eh 

y 

77*00 

Eh 

y 

Eh 

y 

10 

2.963e-02 

— 

7.731e-03 

— 

2.133e-03 

— 

1.347e-03 

— 

20 

9.754e-03 

1.6 

9.713e-04 

3.0 

3.343e-05 

6.0 

2.905e-06 

8.9 

40 

2.452e-03 

2.0 

5.360e-05 

4.2 

3.879e-07 

6.4 

6.511e-09 

oq 

06 

80 

7.076e-04 

1.8 

5.655e-06 

3.2 

1.280e-08 

4.9 

5.687e-ll 

6.8 

160 

1.832e-04 

1.9 

7.544e-07 

2.9 

4.063e-10 

5.0 

2.598e-ll 

1.1 

320 

4.980e-05 

1.9 

9.613e-08 

3.0 

3.017e-ll 

3.8 

4.594e-ll 

-0.8 


Table 2 - and I/°“-errors and observed convergence order v for the M 3 manufactured solution (5.1) with 
optimization gradient tolerance r = 10 “^^. 

using a 100-point Gauss-Lobatto quadrature rule over each spatial cell Ij, and E'^ is approximated by 
taking the maximum over these quadrature nodes. 

The observed convergence order v is defined by 

where for i e {1, 2}, is the error for the numerical solution using cell size Axi, for p S {1, 00}. 

Convergence tables are given in Table 2 for solutions with a tighter optimization tolerance of r = 10“^^. We 
observe that the scheme converges with at least its designed order until the errors are roughly 0{t), where 
errors from the numerical optimization halt the convergence. For many of the solutions on the coarsest 
grids, the convergence is faster than designed, likely because the WENO reconstruction is order 2A: — 1 at 
the cell interfaces. The effects are indeed more pronounced for higher orders. 

In Figure 2 we plot the error of solutions of various orders against their computation time. Here we confirm 
the expectation that for smaller errors, higher-order solutions require less computation time. 

5.2. Plane source 

In this test case we start with an isotropic distribution where the initial mass is concentrated in the middle 
of an infinite domain x S (—00,00): 


tpt=o{x, m) = '(/’floor + < 5 (x) 

where the small parameter ipRoor = 0.5 x 10“® is used to approximate a vacuum. In practice, a bounded 
domain must be used, so we choose a domain large enough that the boundary should have only negligible 
effects on the solution: thus for our final time tf = 1, we take X = [xl,xr] = [—1.2,1.2]. At the boundary 
we set 


V’L(t,M) = V'floor and V'R(i, M) = V’floor 


We set (Tg = 1 and da = 0. 
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Figure 2 - Efficiency of the scheme up to seventh order for the M 3 manufactured solution with optimization 
tolerance r = 10 “^^. 


All solutions here are computed with an even number of cells, so the delta function in the initial condition 
lies on a cell boundary. Therefore we approximate the delta function by splitting it into the cells immediately 
to the left and right. 

Since this is a highly non-smooth problem, the choice of c is not immediately obvious. Indeed, the solution 
contains such strong gradients, so choosing c too small can greatly reduce the accuracy of the simulation. 
After some numerical experimentation, we observed that in the 1000-cell seventh-order simulations the 
value c = 15 generally gave a good trade-off between smoothness in the solution without adding too much 
diffusivity. In the 4000-cell fourth-order simulations, where the non-smoothness is further resolved, c = 10 
gave the best results. 

In Figure 3 and Figure 4, we plot several solutions at the final time for the M3 and M7 models, respectively, 
as well as a reference solution^ of the kinetic equation (2.1). We consider numerical solutions of orders 
fc S {1,4, 7}. The first-order solutions are included to indicate our best guess of the true solutions, and 
they largely agree with those presented in [18] (although we have computed solutions on a much finer grid 
in order to better resolve sharp peaks in the solutions). 

In Figure 3, we present seventh-order M3 solutions. In Figure ??, we see that while the seventh-order 
solution using only the positivity-preserving limiter closely matches the highly-resolved first-order solution, 
some spurious oscillations around x = ±0.25 remain. These oscillations are largely removed by using the 
maximum-principle limiter with c = 15, as we see in Figure ??, though in Figure ?? we see that this solution 
is also somewhat more diffusive. 

The M7 solution to the plane-source problem is more oscillatory. Figure ?? shows that our kinetic scheme 
with fourth-order reconstructions and the maximum-principle limiter performs well. However, seventh-order 
solutions are notably more diffusive, as shown at the leading edge of the solution in Figure ??. 

That the high-order solutions are more diffusive than the true M^v solutions is not necessarily a disadvantage: 


^ The reference solution was computed using a first-order P 199 method on a grid with J = 4000 cells. 
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the more diffusive solutions are actually closer to the reference kinetic solution. Indeed, the equations 
are a spectral method for the original kinetic equation, and thus should not be applied to a non-smooth 
problem like this one without filtering to avoid the Gibbs phenomenon. Exactly how to apply such a hlter 
for moment models kinetic equations is a topic of ongoing research [26, 30], but here it seems that the 
maximum-principle limiter already filters the solution somewhat. 



X 

(a) M 3 with the positivity-preserving limiter 



Figure 3 - The local particle density uq from first- and seventh-order solutions to the plane-source problem 
with the Ms t = 1. First-order solutions were computed with J = 10000 cells and seventh-order with J = 1000 
cells. 


5.3. Source-beam 


Finally we present a discontinuous version of the source-beam problem from [14]. The spatial domain is 
X = [0,3], and 


o-a(a;) 



if a: < 2 

as{x) = 1 

[“ 

if a; < 1 / 

1 ° 

else ’ 

i^O 

if 1 < a; < 2 , 5'(x) = < 

else ^ 


1 

0 


if 1 < a: < 1.5 
else 
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(a) M 7 , fc = 4, c = 10 (b) M 7 , k = 7, c = 15 


Figure 4 - The local particle density uq from first-, fourth- and seventh-order solutions to the plane-source 
problem with the M 7 model at t = 1. First-order solutions were computed with J = 10^ cells; fourth-order with 
J = 4000 cells; and seventh-order with J = 1000 cells. 


with initial and boundary conditions 


= V'floor, 

=/3exp(-7(^ - 1)^) and V'R(i) = V'floor, 
for 7 = 10^, normalization constant /3 = — 1)^))“^, and V’floor = 0.5 x 10“^°. 

M^v solutions for this problem are shown in Figure 5a using J = 150 cells and seventh-order reconstructions, 
using the maximum-principle limiter with c = 1, along with a reference solution.® We see that increasing 
the moment order to iV = 3 qualitatively improves the solution significantly. 

Figure 5b shows some strengths and weaknesses of our scheme for the Mi model, where for this problem the 
shock is the strongest. When compared to a much more finely resolved low-order approximation, the seventh- 
order approximation on a coarse grid nicely fits the features of the solution despite the discontinuous physical 
parameters. However, the incoming beam can only be resolved up to first order (note the piecewise-constant 
approximations behind the beam), resulting in a more diffusive solution at the shock around x = 1. 

The source-beam problem is particularly well-suited to test the oscillation-dampening effects of the limiters. 
In Figure 6 we compare several seventh-order solutions. First in Figure 6a we see in the Mi solution that the 
positivity-preserving limiter does not dampen spurious oscillations while the maximum principle-preserving 
limiter does a better job even for a relatively large value of c. Second, in Figure 6b, we show with the M2 
model that for small values of c, the maximum principle-preserving limiter is too diffusive. 


6. Conclusions and outlook 

In this paper we describe how to implement a kinetic scheme of (in principle) arbitrarily high order for 
entropy-based moment closures of linear kinetic equations in one space dimension. For spatial reconstructions 
we use the well-known WENO method to reconstruct the underlying entropy ansatze using interpolating 


The reference solution was computed using a first-order Pgg method on a grid with J = 2000 cells. 
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(a) Seventh-order Mjv solutions, J = (b) Order comparison for Mi 

150 


Figure 5 - Seventh-order approximation of several Mjv solutions for the sonrce-beam testcase and comparison 
of the seventh-order approximation on a coarse grid with the second-order approximation on a finer grid. 


polynomials, and time integration is performed using multi-step SSP methods. These SSP time integrators 
play a key role in allowing us to give a time-step restriction which guarantees that the moments stay in the 
realizable set. The other key component is a limiter, which not only ensures positivity of the polynomial 
reconstructions on a spatial quadrature set, but also enforces a local maximum principle which dampens 
spurious oscillations in numerical solutions. 

We performed convergence tests with a manufactured solution that included the effects of a space- and 
time-dependent absorption interaction coefficient, and these results validated that the scheme is converging 
at least as fast as expected, and often faster at lower resolutions with higher orders. The convergence tests 
also showed that errors from the numerical optimization routine needed for the angular reconstructions 
limit the overall accuracy of the scheme. This indeed eliminates the benefit of going beyond a certain 
order (depending on the optimization tolerance r). However, our manufactured-solution tests showed that 
increases in efficiency can be obtained before the optimization errors dominate the solution. 

Using the plane-source, an challenging highly non-smooth benchmark problem, we showed that with the 
maximum-principle limiter accurate solutions can be obtained which even limit some of the spurious oscil¬ 
lations due to Gibbs phenomena in the true Mjv solutions, thus pushing our solutions closer to the kinetic 
solution. However, the choice of the parameter c plays an important role in the approximation quality of 
the scheme, and a good value is not always available a priori. With another benchmark problem, the source 
beam, we also demonstrated the benefit of using the maximum principle-preserving limiter. Again, one must 
strike a balance between flattening smooth extrema and dampening spurious oscillations with the choice of 
its parameter. 

Compared to the discontinuous-Galerkin implementation in [4], where realizability preservation is complex 
due to the structure of the set of realizable moments, realizability preservation on the kinetic level is much 
simpler. However, the wider stencils in the WENO reconstruction process will influence the overall paral- 
lelizability of the scheme. 

Future work should continue to work toward practical implementations of entropy-based moment closures. 
Models in two and three spatial dimensions should be implemented, and here a notable challenge is the 
increasing number of angular quadrature points that will be needed. Indeed, our reconstructions are per- 
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1.5 



(a) A closer look at oscillations in the 
Ml solution for different limiter configu¬ 
rations. 


(b) Examining the diffusivity in the 
M 2 model for different limiter configu¬ 
rations. 


Figure 6 - The effects of different limiting in the source-beam test at t = 2.5 for seventh-order solutions with 
J = 150 cells. 


formed at every angular quadrature point, so more efficient WENO techniques will be necessary. Other 
collision models should also be considered. The Laplace-Beltrami operator in the Fokker-Planck equation 
does not fall under the types of collision operators considered here but is an important model for problems 
with forward-peaked scattering. This appears, for example, in important applications such as radiother¬ 
apy. Finally, since we have only considered explicit time-stepping schemes, our time-step restriction scales 
with the mean-free path. Implicit-explicit or asymptotic-preserving schemes should be developed to handle 
moment models near diffusive or fluid regimes without requiring extremely small time steps. 
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